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Abstract.  The  goal  of  this  paper  is  to  present  the  development  of 
a  new  reconstruction  methodology  for  restoring  Magnetic 
Resonance  Images  (MRI)  from  reduced  scans  in  k-space.  The 
proposed  approach  considers  the  combined  use  of  Neural 
Network  models  and  Bayesian  restoration,  in  the  problem  of 
MRI  image  extraction  from  sparsely  sampled  k-space,  following 
several  different  sampling  schemes,  including  spiral  and  radial. 
Effective  solutions  to  this  problem  are  indispensable  especially 
when  dealing  with  MRI  of  dynamic  phenomena  since  then, 
rapid  sampling  in  k-space  is  required.  The  goal  in  such  a  case 
is  to  make  measurement  time  smaller  by  reducing  scanning 
trajectories  as  much  as  possible.  In  this  way,  however, 
underdetermined  equations  are  introduced  and  poor  image 
reconstruction  follows.  It  is  suggested  here  that  significant 
improvements  could  be  achieved,  concerning  quality  of  the 
extracted  image,  by  judiciously  applying  Neural  Network  and 
Bayesian  estimation  methods  to  the  k-space  data.  More 
specifically,  it  is  demonstrated  that  Neural  Network  techniques 
could  construct  efficient  priors  and  introduce  them  in  the 
procedure  of  Bayesian  reconstruction.  These  ANN  Priors  are 
independent  of  specific  image  properties  and  probability 
distributions.  They  are  based  on  training  supervised  Multilayer 
Perceptron  (MLP)  neural  filters  to  estimate  the  missing 
samples  of  complex  k-space  and  thus,  to  improve  k-space 
information  capacity.  Such  a  neural  filter  based  prior  is 
integrated  to  the  maximum  likelihood  procedure  involved  in 
the  Bayesian  reconstruction.  It  is  found  that  the  proposed 
methodology  leads  to  enhanced  image  extraction  results 
favorably  compared  to  the  ones  obtained  by  the  traditional 
Bayesian  MRI  reconstruction  approach  as  well  as  by  the  pure 
MLP  based  reconstruction  approach. 

I.  INTRODUCTION 

A  data  acquisition  process  is  needed  to  form  the  MR  images. 
Such  data  acquisition  occurs  in  the  spatial  frequency  k 
space)  domain,  where  sampling  theory  determines  resolution 
and  field  of  view,  and  it  results  in  the  formation  of  the  k-space 
matrix.  Strategies  for  reducing  image  artifacts  are  often  best 
developed  in  this  domain.  After  obtaining  such  a  k-space 
matrix,  image  reconstruction  involves  fast  multi-dimensional 
Inverse  Fourier  transforms,  often  preceded  by  data 
interpolation  and  re -sampling. 


Sampling  the  k-space  matrix  occurs  along  suitable 
trajectories  [1],  Ideally,  these  trajectories  are  chosen  to 
completely  cover  the  k-space  according  to  the  Nyquist 
sampling  criterion.  The  measurement  time  of  a  single  trajectory 
can  be  made  short.  However,  prior  to  initiating  a  trajectory, 
return  to  thermal  equilibrium  of  the  nuclear  spins  needs  to  be 
awaited.  The  latter  is  governed  by  an  often  slow  natural 
relaxation  process  that  is  beyond  control  of  the  scanner  and 
impedes  fast  scanning.  Therefore,  the  only  way  to  shorten 
scan  time  in  MRI  when  needed,  as  for  instance  in  functional 
MRI,  is  to  reduce  the  overall  waiting  time  by  using  fewer 
trajectories,  which  in  turn  should  individually  cover  more  of  k- 
space  through  added  curvatures.  Although,  however,  such 
trajectory  omissions  achieve  the  primary  goal,  i.e.  more  rapid 
measurements,  they  entail  undersampling  and  violations  of 
the  Nyquist  criterion  thus,  leading  to  concomitant  problems 
for  image  reconstruction. 

The  above  mentioned  rapid  scanning  in  MRI  problem  is 
highly  related  with  two  other  ones.  The  first  is  the  selection 
of  the  optimal  scanning  scheme  in  k-space,  that  is  the  problem 
of  finding  the  shape  of  sampling  trajectories  that  more  fully 
cover  the  k-space  using  fewer  number  of  trajectories.  Mainly 
three  such  alternative  shapes  have  been  considered  in  the 
literature  and  are  used  in  actual  scanners,  namely,  Cartesian, 
radial  and  spiral  [1],  associated  with  different  reconstruction 
techniques.  More  specifically,  the  Cartesian  scheme  uses  the 
inverse  2D  FFT,  while  the  radial  and  spiral  scanning  involve 
the  Projection  Reconstruction,  the  linogram  or  the  SRS-FT 
approaches  [1], 

The  second  one  is  associated  with  image  estimation  from 
fewer  samples  in  k-space,  that  is  the  problem  of  omitting  as 
many  trajectories  as  possible  without  attaining  worse 
reconstruction  results.  The  main  result  of  such  scan 
trajectories  omissions  is  that  we  have  fewer  samples  in  k- 
space  than  needed  for  estimating  all  pixel  intensities  in  image 
space.  Therefore,  there  is  an  infinity  of  MRI  images  satisfying 
the  sparse  k-space  data  and  thus,  the  reconstruction  problem 
becomes  ill-posed.  Additionally,  omissions  usually  cause 
violation  of  the  Nyquist  sampling  condition.  Despite  the  fact 
that  solutions  are  urgently  needed,  in  functional  MRI  for 
instance,  very  few  research  efforts  exist  in  the  literature.  The 
most  obvious  and  simplest  such  method  is  the  so  called  “zero  - 
filling  the  kspace”,  that  is,  all  missing  points  in  kspace 
acquire  complex  values  equal  to  zero.  Subsequently,  image 
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reconstruction  is  achieved  as  usually,  by  applying  the  inverse 
Fourier  transform  to  the  corresponding  k-space  matrix.  Instead 
of  zero -filling  the  k-space  it  might  be  more  advantageous  to 
interpolate  it  by  using  nonlinear  interpolation  procedures,  like 
neural  networks,  as  proposed  by  the  authors  [2],  The 
Bayesian  reconstruction  approach,  developed  by  two  of  the 
authors  [1],  briefly  presented  in  the  next  section  is  another 
alternative  solution.  Both  these  last  two  mentioned  MR 
reconstruction  solutions  yield  good  results  [1,2].  The  main 
contribution,  however,  of  this  paper  is  to  develop  a  novel  MR 
reconstruction  methodology  by  involving  both  Bayesian  and 
Neural  reconstruction  techniques  and  present  its  competence 
and  advantages  over  the  other  rival  approaches. 

n.  THE  BAYESIAN  RECONSTRUCTION  APPROACH. 

The  Bayesian  reconstruction  approach  recently  proposed 
by  two  of  the  authors  [1],  attempts  to  provide  solutions 
through  regularizing  the  problem  by  invoking  general  prior 
knowledge  in  the  context  of  Bayesian  formalism.  The 
algorithm  amounts  to  minimizing  the  following  objective 
function  [1],  by  applying  the  conjugate  gradients  method, 

I  S_-  T II2/  (2a2)  +  (3/2)  £x,y  log  {  a2  +  (xAxy)2  +  (yAxy)2 }  (1) 

with  regards  to  I,  which  is  the  unknown  image  to  be 
reconstructed  that  fits  to  the  sparse  k-space  data  given  in  S. 
The  first  term  comes  from  the  likelihood  term  and  the  second 
one  from  the  prior  knowledge  term  of  the  Bayesian  formulation 
[1],  In  the  above  formula,  T((kx,  ky),(x,y))  =  e~2m(xkx  +  yky) 
represents  the  transformation  from  image  to  kspace  data 
(through  2-D  EFT).  The  second  term  symbols  arise  from  the 
imposed  2D  Lorentzian  prior  knowledge.  xAxy  and  yAxy  are  the 
pixel  intensity  differences  in  the  x-  and  y-  directions 
respectively  and  a  is  a  Lorentz  distribution-width  parameter. 
Assuming  that  P(I)  is  the  prior,  imposing  prior  knowledge 
conditions  for  the  unknown  MRI  image,  then,  the  second  term 
of  ( 1)  comes  as  follows. 

The  starting  point  is  that  P(I)  could  be  obviously  expanded 
into  P(I)=P(Io,o)P(Ii,ol  Io.o)  P(I2,ol  W  Ii,o  )•  •  •  If,  now,  it  is  assumed 
that  the  intensity  Ixy  depends  only  on  its  left  neighbour  ( Ix_ljy 
),  then  the  previous  P(I)  expansion  takes  on  the  form  P(I)  = 
Ilix.y)  P(Ix,yl  Ix-i.y),  provided  that  the  boundaries  are  ignored. 
Next,  we  assume  that  P(Ix,y I  Li,y)  is  a  function  only  of  the 
difference  between  the  corresponding  pixels.  This  difference 
is  written  down  as  xAxy  =  Ixy  -  Ix_1>y.  It  has  been  shown  that  the 
probability  density  function  of  xAxy  is  Lorentzian  shaped  (see 
[1]).  These  assumptions  and  calculations  lead  to  computing 
the  prior  knowledge  in  the  Bayesian  reconstruction  as  in  the 
second  term  of  (1). 

Although  this  Bayesian  reconstruction  approach  tackles 
the  problem  of  handling  missing  samples  in  k-space,  it 
exhibits,  however,  the  disadvantage  that  assumes  the 
existence  of  special  probability  distributions,  given  in  closed 
form  descriptions,  for  representing  the  unknown  ones 


occurred  in  MRI,  which  is  an  issue  under  question.  In  this 
paper  we  attempt  to  remedy  this  problem  by  proposing 
additional  priors  in  the  Bayesian  formulation  in  order  to 
capture  the  probability  distribution  functions  encountered  in 
MRI.  These  priors  are  constructed  through  applying  a 
specifically  designed  Multilayer  Perceptron  (MLP)  neural  filter 
for  interpolating  the  sparsely  sampled  k-space. 

III.  DESIGN  OF  MLP  NEURAL  NETWORK  PRIORS 

The  method  herein  suggested  for  designing  efficient  Priors 
for  the  Bayesian  reconstruction  formalism,  is  based  on  the 
attempt  to  extract  prior  knowledge  from  the  process  of  filling 
in  the  missing  complex  values  in  k-space  from  their 
neighboring  complex  values.  Thus,  instead  of  assuming  a 
Lorentzian  prior  knowledge  to  be  extracted  from  the 
neighboring  pixel  intensities  in  MRI,  as  a  constraint  to  be 
applied  in  the  conjugate  gradient  based  Bayesian 
reconstruction  process,  the  proposed  strategy  doesn’t  make 
any  assumption.  Instead,  it  aims  at  extracting  priors  without 
any  specific  consideration  concerning  the  shape  of  the 
distributions  involved,  by  transforming  the  original 
reconstruction  problem  into  an  interpolation  one  in  the 
complex  domain.  While  linear  interpolators  have  already  been 
used  in  the  literature  [2]  and  nonlinear  estimations  are  well 
established  in  MRI  [2],  ANN  models  offer  several  advantages 
when  applied  as  sparsely  sampled  k-space  interpolators  [2]. 
The  methodology  to  extract  prior  knowledge  by  applying  the 
ANN  filters  in  MRI  reconstruction  is  described  in  the 
following  paragraphs. 

Stepl.  We  compile  a  set  of  R  representative  N  X  N  MRI 
images  with  k-space  matrices  completely  known,  which 
comprise  the  training  set  of  the  MLP  interpolators. 
Subsequently,  we  scan  these  matrices  following  the  specific 
sampling  schemes  mentioned  above  and  then,  by  randomly 
omitting  trajectories  the  sparse  k-spaces  are  produced. 

Step2.  The  original  k-space  matrix  as  well  as  its  corresponding 
sparse  k-space  matrix  associated  with  one  N  X  N  MRI  training 
image,  is  raster  scanned  by  a  (2M+1)  X  (2M+1)  sliding 
window  containing  the  associated  complex  kspace  values. 
The  estimation  of  the  complex  number  in  the  center  of  this 
window  from  the  rest  of  the  complex  numbers  comprising  it  is 
the  goal  of  the  proposed  interpolation  procedure.  Each 
position  of  this  sliding  window  is,  therefore,  associated  with  a 
desired  output  pattern  comprised  of  the  complex  number  in 
the  original  k-space  corresponding  to  the  window  position, 
and  an  input  pattern  comprised  of  the  complex  numbers  in  k- 
space  corresponding  to  the  rest  (2M+1)  X  (2M+1)  -1  window 
points. 

Step3.  Each  such  pattern  is  then,  normalized  according  to  the 
following  procedure.  First,  the  absolute  values  of  the  complex 
numbers  in  the  input  pattern  are  calculated  and  then,  their 
average  absolute  value  lzaverl  is  used  to  normalize  all  the 
complex  numbers  belonging  both  in  the  input  and  the  desired 
output  patterns.  That  is,  if  ?  is  such  a  number  then  this 


normalization  procedure  transforms  it  into  the  Z|/lzaverl.  In  the 
case  of  test  patterns  we  apply  the  same  procedure.  That  is,  the 
average  absolute  value  lzaverl  for  the  complex  numbers  z,  of  the 
test  input  pattern  is  first  calculated.  Then,  the  normalized 
complex  values  z/lzaverl  feed  the  MLP  interpolation  filter  to 
predict  the  sliding  window  central  normalized  complex  number 
znormCentre-  The  corresponding  unnormalized  complex  number  is 

Simply  Z  centre  *  lz^vej.1. 

Step4.  The  next  step  is  the  production  of  training  patterns  for 
the  MLP  interpolators  and  their  training  procedure.  To  this 
end,  by  randomly  selecting  sliding  windows  from  the 
associated  k-spaces  of  the  R  training  images  and  producing 
the  corresponding  input  and  desired  output  training  pairs  of 
patterns,  as  previously  defined,  we  construct  the  set  of 
training  patterns.  The  assumption  underlying  such  an 
approach  of  training  MLP  interpolators  is  that  there  are 
regularities  in  every  k-space  sliding  window,  the  same  for  any 
MRI  image,  to  be  captured  by  the  MLPs  without  any  prior 
assumption  for  the  probability  distributions.  MLP  training 
follows  by  applying  the  conjugate  gradient  technique  of 
Polak-Ribiere. 

Step5.  After  training  phase  completion,  the  MLP  filter  has 
been  designed  and  can  be  applied  to  any  similar  test  MRI 
image  as  follows.  To  this  end,  the  (2M+1)  X  (2M+1)  sliding 
window  raster  scans  the  sparse  k-space  matrix  associated  with 
this  test  image,  starting  from  the  center.  Its  central  point 
position  moves  along  the  perimeter  of  rectangles  covering 
completely  the  k-space,  having  as  center  of  gravity  the  center 
of  the  k-space  array  and  having  distance  from  their  two 
adjacent  ones  of  1  pixel.  It  can  move  clockwise  or 
counterclockwise  or  in  both  directions.  For  every  position  of 
the  sliding  window,  the  corresponding  input  pattern  of 
(2M+1)  X  (2M+1)  -  1  complex  numbers  is  derived  following 
the  above  described  normalization  procedure.  Subsequently, 
this  normalized  pattern  feeds  the  MLP  interpolator.  The 
wanted  complex  number  corresponding  to  the  sliding  window 
center,  is  found  as  Centre  =  AiLP°ut  *  lzaverl,  where  Zmlp01"  is  the 
MLP  output  and  lzaverl  the  average  absolute  value  of  the 
complex  numbers  comprising  the  unnormalized  input  pattern. 
For  each  rectangle  covering  the  k-space,  the  previously 
defined  filling  in  process  takes  place  so  that  it  completely 
covers  its  perimeter,  only  once,  in  both  clockwise  and 
counterclockwise  directions.  The  final  missing  complex  values 
are  estimated  as  the  average  of  their  clockwise  and  counter¬ 
clockwise  obtained  counterparts.  The  outcome  of  the  MLP 
filter  application  is  the  reconstructed  test  image,  herein  named 
MLP_Img  (equation  (2)).  Its  difference  from  the  image  I<u 
obtained  during  the  previous  step  of  conjugate  gradient 
optimization  in  the  Bayesian  reconstruction  formula  (1), 
provides  the  neural  prior  to  be  added  for  the  current 
optimization  step. 


IV.  INCORPORATION  OF  NEURAL  NETWORK  PRIOR 
KNOWLEDGE  INTO  THE  BAYESIAN  PRIOR 


Following  the  5  steps  above  we  can  formulate  the 
incorporation  of  MLP  priors  to  the  Bayesian  restoration 
process  as  follows. 

•  Design  the  MLP  Neural  filter  as  previously  defined 

•  Consider  the  Bayesian  reconstruction  formula  (1).  The 
image  to  be  optimized  is  I  given  the  k-space  S.  The  initial 
image  in  the  process  of  conjugate  gradient  optimization  is  the 
zero-filled  image.  At  each  step  t  of  the  process  a  different  I(t) 
(the  image  at  the  t  step,  that  is,  the  design  variables  of  the 
problem)  is  the  result.  Based  on  fig.  1,  by  applying  the  MLP 
filter  on  the  original  sparse  kspace,  but  with  the  missing 
points  initially  filled  by  the  FFT  of  I(t)  (in  order  to  derive  the 
I(t)  k-space)-  and  afterwards  refined  by  the  MLP  predictions, 
we  could  obtain  the  difference  I(t)  -MLP_Img(t)  as  the  Neural 
Prior. 

•  The  Neural  Network  (NN)  Prior  form,  therefore,  is: 

^1  MLP _lmg{t\x,y)  - I{t)(x ,  y )  I  (2) 

x,y 

where,  MLP_Imgm(x,y)  is  the  NN  estimated  pixel  intensity 
in  image  space  (NN  reconstructed  image:  Inverse  FFT  of  NN 
completed  k-space)  at  step  t  and  I^fey)  is  the  image  obtained 
at  step  t  of  the  conjugate  gradient  optimization  process  in  the 
Bayesian  reconstruction 

•  The  proposed  Prior  in  the  Bayesian  reconstruction  is 
given  as 

Final  Prior  =  Lorentzian  Bayesian  Prior  +  a*  NN_Prior  (3) 

•  That  is,  the  optimization  process  f”  is  attempted  to  be 
guided  by  the  MLP_Img  ®  produced  by  the  NN 
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Figure  1.  The  difference  between  the  image  to  be  optimized  I(t)  and  the 
MLP  reconstructed  image  MLP_Img(t)  constitutes  the  neural  prior. 


V.  EVALUATION  STUDY  AND  CONCLUSIONS 

An  extensive  experimental  study  has  been  conducted  in  order 
to  evaluate  the  above  defined  novel  Bayesian  reconstruction 
methodology.  All  the  methods  involved  have  been  applied  to 
a  large  MRI  image  database  downloaded  from  the  Internet, 
namely,  the  Whole  Brain  Atlas  http://www.med.harvard.edu/ 
AANLIB/home.html  (copyright  ©  1995-1999  Keith  A. 
Johnson  and  J.  Alex  Becker).  We  have  used  10  images 
randomly  selected  out  of  this  collection  for  training  the  MLP 
filters,  and  10  images,  again  randomly  selected  for  testing  the 
proposed  and  the  rival  reconstruction  methodologies.  All 
images  are  256  by  256.  Their  kspace  matrices  have  been 
produced  applying  the  2D  FFT  to  them.  Radial  trajectories 
have  been  used  to  scan  the  resulted  256  X  256  complex  k- 
space  arrays.  4  X  256  =  1024  radial  trajectories  are  needed  to 
completely  cover  such  k-spaces.  In  order  to  apply  the 
reconstruction  techniques  involved  in  this  study,  each  k 
space  has  been  sparsely  sampled  using  128  only  radial 
trajectories.  Regarding  the  sliding  window  raster  scanning  the 
k-space,  a  5  X  5  window  was  the  best  selection. 

Concerning  the  MLP  filter  architecture,  the  48-10-2  one  was 
found  to  be  the  best.  This  MLP  has  been  trained  using  3600 
training  patterns.  The  compared  reconstruction  techniques 
involved  in  this  study  are:  the  proposed  novel  Bayesian 
reconstruction  approach,  the  traditional  Bayesian 
reconstruction  technique  as  well  as  the  MLP  filtering 
interpolation  technique.  Moreover,  the  simplest 
“interpolation’'  approach,  namely  filling  in  the  missing 
samples  in  k-space  with  zeroes  and  then,  reconstructing  the 
image,  has  been  invoked.  All  these  methods  have  been 
implemented  using  the  MATLAB  programming  platform. 

Concerning  the  measures  involved  to  quantitatively 
compare  reconstruction  performance,  we  have  employed  the 
usually  used  Sum  of  Squared  Errors  (SSE)  between  the 
original  MRI  image  pixel  intensities  and  the  corresponding 
pixel  intensities  of  the  reconstructed  image.  Additionally, 
another  quantitative  measure  has  been  used,  which  expresses 
performance  differences  in  terms  of  the  RMS  error  in  dB  [2]: 
lambda=(image_recon(:)'*image_orig(:))/(image_recon(:)'*im 
age_recon(:));residu=image_orig-lambda*image_recon; 
dB=10*logl0((image_orig(:)'*image_orig(:))/(residu(:)'* 
residu(:))); 

The  quantitative  results  obtained  by  the  different 
reconstruction  methods  involved  are  outlined  in  table  1 
(average  SSE  and  RMS  errors  for  the  10  test  MRI  images). 
Concerning  reconstruction  performance  qualitative  results,  a 
sample  is  shown  in  figure  2.  Both  quantitative  and  qualitative 
results  clearly  demonstrate  the  superiority  of  the  proposed 
Bayesian  reconstruction  methodology  embedding  MLP 
filtering  based  prior  knowledge,  in  terms  of  MRI  image 
reconstruction  performance  over  the  other  three  rival 
methodologies  (simple  Bayesian,  MLP  filter  and  zero -filled 
reconstructions).  Future  trends  of  our  research  effort  include 


implementation  of  the  3-D  Bayesian  reconstruction  with  NN- 
priors  for  f-MRI  as  well  as  applications  in  MRI  image 
segmentation  for  tumor  detection. _ 


Method 

SSE  (average  in 
the  10  test  MRI 
images) 

dB  (average  in  the 
10  test  MRI 

images) 

Proposed  Bayesian 
reconstruction  with 

NN  Prior 

2.85  E3 

16.67 

Bayesian 

reconstruction 

3.40e3 

15.92 

MLP  filtering 

3.30E3 

15.98 

Zero-filling 

reconstruct. 

3.71E3 

15.26 

Table  1.  The  quantitative  results  with  regards  to 
reconstruction  performance  of  the  various  methodologies 
involved 


Fig.  2  From  left  to  right  and  top  to  down:  The  proposed 
Bayesian  with  NN  prior,  the  sparsely  sampled  k-space 
(nr=128)-zerofilled  image  reconstruction,  the  MLP  filtering  and 
the  traditional  Bayesian  reconstruction  results.  The  Test 
Image  illustrates  a  brain  slice  with  Alzheimer’s  disease 
(http://www.med.harvard.edu/  AANLIB/cases/case3/mrl- 
tcl/020.html). 

References 

[1]  G.H.L.A.  Stijnman,  et  al.:  MR  Image  Estimation  from 
Sparsely  Sampled  Radial  Scans  "Proc.  ProRISC/IEEE  Benelux 
on  Circuits,  Systems  and  Signal  Processing,  Mierlo  (The 
Netherlands),  1997" ,  603-61 1 

[2]  M.  Reczko,  D.A.  Karras,  et  al.,  “Neural  Networks  in  MR 
Image  Estimation  from  Sparsely  Sampled  Scans”,  MLDM’99, 
Leipzig,  1999,  in  Lecture  Notes  In  Artificial  Intelligence  1715, 
Springer,  pp.  75-86, 1999. 


